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This  wjrk  was  motivated  during  a study  of  radiative  transfer 
through  aerosols  when  it  was  realized  that  spherical  aerosol  cress 
sections  and  phase  functions  only  could  be  provided  by  Mie  theory. 
Because  many  aerosols  are  irregular  in  shape,  another  formalism  is 
needed  to  provide  the  optical  parameters  needed  for  radiative  transfer 
calculations.  The  literature  of  previous  efforts  on  the  scattering  of 
electromagnetic  waves  from  irregular  particles  is  extensive  , although 
many  of  these  calculations  were  directed  either  to  various  approximate 
approaches  or  to  the  problem  of  conducting  materials  only.  An  example 
of  both  these  limitations  is  the  application  of  least-square^  techni- 
ques to  match  electric  fields  at  the  surface  of  a conductor.'  How  the 
least-squares  solution  is  related  £o  the  boundary  conditions  is  unclear. 

Several  recent  efforts  using  integral  equation  formulations 
deserve  special  attention  since  not  only  are  the  boundary  conditions 
properly  accounted  for,  but  numerical  results  are  presented.  Referen- 
ces (h)  and  (5)  are  based  on  the  scalar  rather  than  the  vector  wave 
equation.  Reference  (6)  is  completely  general,  but  was  only  applied 
to  conducting  particles.  Another  entirely  different  effort  of  note„ 
which  will  be  addressed  later  was  that  of  Chylek,  Grams  and  Pinnick 
(hereafter  Chylek  et  al. ) Chylek  et  al.  modify  the  Mie  theory  to 
account  for  particle  irregularity  on  the  empirical  grounds  that  the 
phase  functions  of  irregular  particles  do  not  exhibit  glories.  As  a 
result,  they  modify  the  resonance  behavior  of  the  field  expansion 
coefficients  to  model  irregular  particle  behavior.  Chylek  et  al.  have 
presented  calculations  that  enjoyed  better  agreement  with  experimental 
data  than  did  Mie  theory. 
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In  this  paper,  the  scattering  of  a plane  electromagnetic 
wave  by  an  irregular  but  cylindrically  symmetric  dielectric  particle 
is  solved  by  generalizing  the  original  work  of  Mie.  ’ As  in  the 
original  work,  the  expansion  coefficients  of  the  scattered  and  inter- 
ior (to  the  dielectric  particle)  electric  and  magnetic  fields  ore 
obtained  by  satisfying  the  boundary  conditions.  Because  of  the  par- 
ticle irregularity,  the  expansion  coefficients  are  coupled  and  form  a 
set  of  linear  equations  that  must  be  solved  numerically.  Both  the 
integral  and  differential  equation  approaches  are  equivalent  in  analogy 
to  the  calculation  of  the  scattering  amplitude  in  quantum  mechanics 
by  solving  either  the  Schrodinger  or  the  Lippman-Schwinger  equation." 
Further,  both  approaches  develop  into  solutions  of  linear  algebraic 
equations.^- 

It  must  be  emphasized  that  both  References  ( U ) — ( 6 ) and  this 
paper  are  limited  in  applicability  to  "slightly"  irregular  particles 
whose  deviation  from  a smooth  geometric  shape  is  much  less  than  the 
wavelength.  This  limitation  on  solving  the  scattering^problern  for 
"arbitrary"  irregular  particles  is  reviewed  by  Millar.^  In  essence 
the  limitation  arises  from  the  fact  that  the  scattered  field  solutions 
incorporate  radial  functions  that  become  outgoing  spherical  wave 
representations  in  the  far  field,  but  not  functions  that  represent 
scattering  of  light  from  one  area  on  the  surface  of  the  particle  to 
another  area.  Consequently,  the  problem  of  multiple  scattering 
between  areas  on  the  particle  have  to  be  neglected.  While  the  over- 
all problem  may  be  handled  by  dividing  the  particle  into  subparti cle^ 
and  satisfying  the  boundary  conditions  for  and  between  all  particles, 
this  paper  will  be  limited  to  only  "slightly"  irregular  particles. 

The  differential  equation  approach  is,  of  course,  also  valid 
for  "slightly"  irregular  particles  without  cylindrical  symmetry,  but 
these  calculations  require  more  computer  core  space  than  is  commonly 
available.  In  Section  II  of  this  paper,  the  Mie  solution  for  a 
spherical  particle  is  reviewed  by  developing  the  differential  equation 
formalism.  Section  III  is  an  extension  of  this  formalism  to 
cylindrically  symmetric  irregular  particles  by  developing  the  four 
boundary  condition  equations  into  UL  equations  that  may  be  solved 
numerically  for  the  expansion  coefficients  to  arbitrary  accuracy.  In 
Section  IV,  comparisons  yi,th  exact  calculations  for  spheres  and 
ellipsoids  of  revolution  ^ are  presented  and  the  computation  problems 
for  more  general  particles  are  discussed.  The  effects  of  irregularity 
on  the  scattering  phase  functions  for  single  particles  are  reviewed, 
and  calculations  performed  for  comparison  with  Reference  (7)  are 
presented  and  discussed. 

II.  PRELIMINARY 

The  Hertz  vectors  II  and  E are  related  to  the  electric  and 
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magnetic  fields  E and  H by 


2_2. 


E = -tkVxn  + k in  E + V(V*E) 


and 


H = k?m2n  + V(V«n  ) - -tkm2VxE 


(1) 

(2) 


where  m , the  square  of  the  complex  refractive  index  is 


m = e + l»it*.o/u) 


(3) 


e is  the  dielectric  constant,  a is  the  conductivity,  u>  is  the  radial 
frequency,  k is  the  wave  number  of  the  monochromatic  incident 
light,  and  u>  = ck  (c  = speed  of  light. ) 

^ Hie  requirement  that  E and  H be  transverse  as  r-*“>  (the  far 
field)',  implies  that  only  the  radial 'components  of  n and  E need  be 
considered. 


n = fir 
E = fa 


(M 

(5) 


where  i*  is  the  radial  unit  vector,  and  u and  a are  known  as  Debye 
potentials.  The  vector  wave  equation  satisfied  by  n and  a reduces 
to  a scalar  equation 


(V2+  k2m2) 


= 0 


(6) 


Equations  (l),  (2),  (h),  and  (5)  may  be  used  to  write  the 
components  of  E and  H in  terms  of  the  Debye  potentials  as 
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If  the  incident  light  is  u pliute  wave  propagating  along  the  z axis, 
then 

Ex  = expU'kz)  (13) 


Hy  = texp(tkz)  (1*0 

which  give  rise  to  Debye  potentials  of  the  fora 

. -2c  -i-\  2f+l  . v „1,  „JcosU)l  ,,,v 

n1)  = k l 1 (151 T ^Ur)  VCOSe)[SinU)J  {U) 

where  the  superscript  "i"  indicates  the  incident  wave.  The  general 

solution  of  Equation  (6)  has  the  form 
L,M 

S = l {a^^(kmr)+3^n^(krnr)}P^(cos0){y^in(miji)+6mcos(m<|>)}(l6) 

whereip^(x)  = xj^Cx),  h|.(^)=  x/^p(x)^  /»(x)  and  </^(x)  are  spherical 
Bessel  functions  iuid  is  the  m 1 Legendre  polynomial  of  order  i . 
•l'he  interior  ("w")  and  scattered  ("s")  potentials  may  be  written  from 
Kq.  (16)  by  noting  from  Eq.  (15)  that  only  m=l  terms  will  appear. 
Further,  the  interior  potentials  must  converge  to  zero  as  r-*-0  so 
that  =0,  and  the  scattered  potentials  must  converge  to  spherical 
waver,  as  r-*»  , so  that  3/.=  ta^.  '^'his  latter  condition  is  equivalent 
to  replacing  ^ iuid  n with  i}£=  x/l^  » where  /t^is  the  spherical  Hankel 
function  of  the  first  kind. ‘ These  potentials  may  be  written  as 

£)■  p>“0)  (17) 

!•)■  '"llipS]  >’f(kr)  rf(cos8)  (18) 

b 

For  a sphere,  the  transverse  continuity  boundary  condition 
requires  that  E , E H^,  and  H be  continuous  across  the  boundary 
of  the  sphere.  Inspection  of  Eqs.  (Y)-(ll’)  reveals  that  this 
boundary  condition  will  be  satisfied  if 


3 / i s > i 3 w | 

3r  'r=u  3r  1 r=a 


3 / i 8 . | 3 _w i 

3r'  ’r=a  3r  1 r=a 


FOWLER 


/ i 8 v | 2 Vi 

(o  +0  )|  = in  a I 


r=a 


(h^+ttS)  | = nW| 

1 r=a  r=e 


(21) 

(22) 


where  a is  the  radius  of  the  sphere.  Equations  (l9)-(22)  may  be 
reduced  by  using  the  orthogonality  of  the  , cos  ($),  and  c in ( <J> ) 
and  the  explicit  forms  of  the  Debye  potentials  to 


/ 1 ? P+i  i 1 

1 KUI y^(*a)  + ^’(ka)  = ;aA'(ta)  (23) 

£ffe^£(ka)  + vl'(ka)  = W(kma)  (2U) 

^-1  zH+l)^£rka^  + c£p^ka)  = a£^(kma)  (25) 

/_1  + d^p^(ka)  = ib^^(kma)  (26) 


where  i|i»'(x)  =3^^(x)/3x.  Equations  (23)-(26)  are  commonly  solved  for 
c^  and  d^  and  are  called  the  Mie  solution.  These  c«  and  d^  are  then 
used  to  calculate  the  cross  sections  and  phase  functions. 

III.  CALCULATIONS 

The  radius  of  an  irregular,  but  cylindrically  symmetric 
particle  may  be  defined  by 

r = Ys  P (cosO)ii  a(0).  (27) 

••  n n 

n-v 

The  normal  vector  to  the  surface  is 
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The  boundary  conditions  are 
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ftXA”>r=u(0)  = 0 (3°) 

where 

AE  = E1  + ES  - EW 
AH  = H1  + HS  - HW. 

These  boundary  conditions  may  be  explicitly  written  as 


AV 

0 

(33) 

AH  = 
<?> 

0 

(3M 

AE0- 

AE  = 0 

30  r 

(35) 

AH0- 

31na( 0 ) , _ 

30  AHr  " 0> 

(36) 

(31) 

(32) 


the  AE  and  AH  term 
r r 

spherical  nature  of  the  particle. 


Equations  (33)— (36)  differ  from  the  Mie  boundary  conditions  only  in 

in  Eqs.  (35)  and  (36)  that  arise  from  the  non 
These  equations  cannot  be  simpli- 
fied in  the  same  manner  as  the  Mie  solution  by  exploiting  the  ortho- 
gonality of  the  angular  functions  because  the  radius  of  the  particle 
is  itself  a function  of  angle  0.  Equations  (31)  and  (32)  may  be 
used  to  rewrite  Eqs.  ( 33 ) — ( 36 ) as 


„s  „w 
4>  “ > ~ 


-E 


H 


s 


* 


W 1 

h:  = -h: 


4 - < - * 


i 3lna( 0)  i 
-E0  30  Lr 


(37) 

(38) 

(39) 


s w 31na(0) . s w.  _ 
H0  “ H0  " DO  (Hr~V  ‘ 


i 31.na(.0l  i 
0 30  r 


(»*0) 


The  individual  components  of  the  electric  and  magnetic  fields  may  be 
calculated  in  terms  of  the  Debye  potentials  by  use  of  Eqs.  (7)-(l2) 
and  substituted  into  Eqs.  (37)-(**0).  These  resulting  equations  are 
functions  only  of  the  angles  0 and  incorporating  ^1-  unknowns;  a^,, 
b»,  c»,..and  d~.  By  calculating  the  overlap  of  Eqs.  (37)  and  (39)  1 
with  P»,(cos0)  eos(if>)  and  Eqs.  (38)  and  (In'))  with  Pp(cos0)  sin($), 
UL  equations  in  unknowns  may  be  fonned.  It  may  be  noted  that  the 
$ integrations  may  be  performed  immediately . The  resulting  U L 
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equations  may  be  expressed  in  figurative  form  as 
L L 

I Bf  Pi  GP  = I FP»Pi 
1=1  4 C /=! 


(Ul) 


where 
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Each  element  A»»  , is  an  L by  L array  whose  entries  are  the  integrals 
given  for  n = I,  6?11  (by  steps  of  five)  by 


where 


A..,  = [ f (a(  0))  f~-na^-  P*(  COS0  )P*  (cose)  dcose  (1»5) 

££  »n  J n -t  4. 

ae+A-U  ( 

£(£+l)  £ 


f,  =/_1  T^^da)  ++."(ka)} 


7 


Bowler 


and 


where 


= (^(kma)+  ♦jj(tana)} 

(U6) 
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«11  = pl  {ka)- 


Equation  ( Ul ) may  be  solved  numerically  Cor  the  a^,  b^,  c»,  and  d^, 
and  the  cross  sections  and  phase  functions  may  be  calculated  from  the 
c^  Rnd  d^  using  the  same  formulas  as  before. 


IV.  CALCULATIONS 

A computer  code  was  generated  to  perform  the  calculations 
described  in  the  previous  section.  The  integrals  A„ „ , ..  Eqs.(U5)- 

(51)  , were  performed  using  a Gauss-Legendre  Quadrature?  Equation 
(Ul)  was  solved  using  j Gjjiss-Jordan  elimination  routine  modified  to 
treat  complex  entries.  ’ The  cro^s  factions  and  phase  functions 
were  calculated  in  the  usual  manner.  ’ 

The  theoretical  development  advanced  in  Section  II  is,  of 
course,  also  valid  for  particles  not  cylindrica-lly  symmetric  with 
respect  to  the  incident  plane  wave.  At  this  time,  however,  an 
extension  of  the  theory  and  the  computer  code  to  consider  these  more 
general  geometries  is  not  feasible  due  to  computer  limitations.  If  a 
general  particle  or  incident  direction  were  to  be  considered,  Eqs.  (l?) 
and  (l8)  would  include  terms  cos(rmj))  and  sin(m<{i)  ,m=  O^to  L.  gTne 
size  of  the  array  B would  increase  to  approximately  Ul  by  Ul  since 
each  current  entry  would  increase  by^terms  m and  m'.  Arrays  G and  F 
would  increase  in  size  from  UL  tg  Ul  . Thus  jbhe  total  storage 
required  would  increase  from  16L  + 8L  to  l6L  + 8l2.  For  La»L2»L, 


this  is  iui  approximate  increase  by  a factor  of  l?.  The  increase  in 
computer  memory  required  for  operational  code  (as  compared  to  array 
storage)  should  be  approximately  proportional.  The  present  code 
requires  2 x 8*  words  to  initiate  compilation,  so  that  L-  ~ 8-i  words 
and  a general  code  will  require  approximately  2 x &l>  words  to  initiate 
compilation.  Computer  memory  of  this  extent  is  available  only  on  a 
very  few  machines  so  that  the  implementation  is  not  feasible  at  this 
time. 

It  should  be  noted  that  while  this  code  requires  a large 
amount  of  computer  memory  for  initial  compilation  (.'?  x 8-’  words,)  this 
does  not  restrict  the  utility  of  the  code  as  much  as  might  be  expected . 
This  code  is  normally  exercised  on  the  MICOM  CPC  6600  computer  using 
the  SCOPE  3.1*. 2 compiler.  This  is  a two  pass  compiler  that  performs 
some  optimization  of  the  compiled  code.  In  this  case,  the  optimization 
is  significant  (is  the  computer  memory  required  for  the  compiled  code 
is  only  6.5  x 8‘*  words.  Additionally,  execution  of  the  code  is 
accomplished  in  about  15  seconds  for  a nominal  particle  of  up  to  Id 


terms  in  the  radius  (Eq.(lV)  and  (l8))  and  Mie  parameter  of  value 


although  compilation  of  the  code  requires  about  60  seconds. 


i bus 


line , 
the 


compiled  code  may  be  used  to  perform  calculations  with  only  moderate 
demands  on  computer  memory  and  execution  time  in  cont rust  to  recom- 
pilation of  the  code  for  each  particle  to  be  considered. 

The  code  was  checked  by  exercising  it  for  various  sizes 
of  spherical  particles  and  various  complex  refractive  indices.  These 
results  ( Cn , d,. .cross  sections  and  phase  functions)  were  compared  with 
the  equivalent  quantities  calculated  using  a standard  Mie  cede.'1  The 
c.  and  d for  all  calculations  agreed  to  within  the  round-off  error 
or  the  machine  (fifteenth  figure.) 

An  additional  check  was  made  by  making  calculations  for 
proluto  spheroids  of  semimajor  to  semiminor  axis  ratio  of  two.  These 
calculations  were  compared  with  exact  calculations.  *-  There  were 
minor  differences  between  the  approximate  calculations  of  this  paper 
and  the  exact  calculations  in  terms  of  the  fine  structure  of  the 
phase  functions,  but  the  positions  of  relative  maxima  and  minima  and 
the  relative  magnitude  of  the  curves  agreed  very  well.  These  calculat- 
ed phase  functions  are  shown  in  Figure  (l).  While  a precise  estimate 
of  accuracy  is  difficult,  a comparison  of  calculated  values  with 
values  fran  the  exact,  curves  at  10°  increments  gave  an  error  bound  of 


?%  for  value  of  c = 2na/.75A=  1,3, 5.7  (a  = semimajor  axis)  with  an 
average  error  of  approximately  8 x 10“**  %.  Because  some  of  this 
error  may  be  attributed  to  computer  noise  and  the  truncation  error  of 
the  Gauss-Legendre  integration,  this  error  bound  seems  to  be  reason- 
able. The  agreement  shows  that  the  code  is  use  Ail  despite  the 
assumption  of  the  validity  of  the  Rayleigh  Hypothesis. 

The  modification  to  Mie  theory  advanced  by  Chylek  et  al. 
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wus  examined.  In  this  modification,  based  on  the  consideration  that 
glories  are  not  observed  in  experimental  data,  the  particle  irregular- 
ity is  modeled  by  truncation  of  the  resonances  in  the  c p and  d^  in  the 
region  2nr/X-£when  l> 3 in  a standard  Mie  code.  While  this  modification 
is  attractive  because  it  offers  a simple,  efficient  computational 
treatment  of  non  spherical  particle  scattering,  the  validity  of  the 
approximation  has  been  questioned.  0 

As  a starting  point  in  considering  the  Chylek  et  al.  modi- 
fication, the  cross  sections  and  phase  functions  for  several  variations 
of  particle  size  and  shape  were  calculated.  In  terms  of  the  Mie  para- 
meter, 

Xq  = 2irr/X , (52) 

particles  of  shape  function 

x(G)  = xQ  {1+  .1  Pn(Q)>  (53) 

for  x = 1,3, 5, 7,9  and  n = 2,3,4  were  considered.  A refractive  index 
of  1.5  + 0-t  was  used  in  this  and  all  subsequent  calculations.  Cross 
sections  both  greater  than  and  less  than  the  equivalent  Mie  cross 
sections  for  the  same  value  of  x0  were  observed.  The  forward  scatter 
and  backscatter  varied  similarly,  in  some  cases  being  greater  and  in 
other  less.  The  structure  of  the  phase  functions  were  similar  but 
not  identical  to  the  Mie  phase  functions.  Although  it  is  difficult 
to  draw  concise  trends,  it  was  noted  that  in  general  differences 
increased  with  increasing  n and  decreased  with  increasing  xQ. 

Additionally,  the  effect  of  varying  the  strength  of  the  de- 
formation was  examined.  Shape  functions  of  the  form 

x(e)  = xQ  {i  + anPn(e)}  (54) 

for  n = 2 and  3 , an  = .05,  .1,  and  .15  were  considered.  Two  effects 
of  this  variation  of  an  were  noted;  the  values  of  the  relative  maxima 
and  minima  varied  by  as  much  as  a factor  of  five  although  the  positions 
of  the  maxima  and  minima  varied  only  slightly,  and  the  amount  of  near 
forward  scattering  ( scattering  angle  less  than  30°)  varied  by  as  much 
as  a factor  of  five  about  the  equivalent  Mie  curve.  This  latter 
effect  would  seem  significant  in  terms  of  the  experimental  data 
reported  by  Chylek  et  al.  The  significance  of  the  other  effects  is 
more  difficult  to  assess  due  to  the  averaging  process  performed  when 
the  polydisperse  phase  functions  are  calculated. 

The  behavior  of  the  c^  and  d^  in  the  region  x0=£  was  studied 
to  examine  the  validity  of  the  modification  of  Chylek  et  al.  Shape 
function  of  the  form  of  Eq.  (53)  for  xQ  = 3 to  4.5  by  steps  of  .1,  and 
n = 2,3,4, 5 were  considered.  The  real  and  imaginary  parts  of  C3  and 
were  compared  with  the  C3  and  d3  calculated  using  Mie  theory.  In 
seven  of  these  eight  calculations,  the  resonance  curve  shape  is 
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preserved  although  it  narrows  by  as  much  as  a factor  of  two.  Addi- 
tionally, the  center  of  the  resonance  curve  shifts,  apparently  without 
trend,  although  the  resonance  curves  (real  parts)  for  the  irregular 
particles  always  fall  under  the  Mie  resonance  curves.  In  one  case, 
that  of  x(0)  = xb{l+.lP3( 0) } , the  resonance  curve  is  replaced  by  a 
slowly  varying  curve  contained  under  the  Mie  curve  (rela  part)  whose 
real  parts  have  a mean  of  approximately  .3  and  whose  imaginary  parts 
have  a mean  of  approximately  zero.  In  all  other  cases,  the  resonance 
curves  narrow  and  fall  under  the  Mie  curves.  The  extent  of  narrowing 
appears  to  be  approximately  proportional  to  n.  On  the  basis  of  these 
calculations,  the  modification  of  the  c^  and  d^  in  the  resonance 
region  is  not  substantiated,  but  the  numerical  results  of  Reference 
(7)  might  be  justified  in  view  of  the  narrowing  of  the  resonance 
curves. 

Calculations  were  also  performed  for  one  of  the  polydisperse 
aerosol  distributions  reported  by  Chylek  et  al.  (i.e.  the  first  KC1.) 
Nine  phase  functions  using 


x( 0)  = x0(I  + .05  P2(0)) 

x(0)  = xc{l  + .1  P2(0)) 

x( 0)  = x0{l  + .15  P2(0)) 

x(0)  = Xq{ 1 - .1  P2(0)  + .1  P^O)} 

were  calculated  and  combined  linearly  with  adjustable  coefficients 
calculated  using  regression  techniques  to  far  the  data.  Two  different 
fitting  techniques  were  used;  one  using  all  the  data,  and  another 
using  only  those  data  for  scattering  angles  greater  than  sixty  degrees. 
A representative  fit  for  each  scheme  is  shown  in  Figure  (2).  These 
fits  appear  to  be  as  good  as  that  reported  by  Chylek  et  al.  It  is 
noted  that  the  fit  is  not  sensative  to  the  choice  of  combinations  of 
phase  functions  used. 

The  insensitivity  of  the  fits  to  variation  in  the  phase 
functions  used  in  the  fitting  calculations  seems  to  indicate  that  a 
more  elaborate  calculation  for  more  irregular  shaped  particles  averaged 
over  orientations  would  not  significantly  improve  the  results  in  terms 
of  the  cost  of  such  a calculation.  In  the  same  vein,  the  Chylek  et  al. 
modification  seems  to  be  adequate  for  many  calculations.  It  appears 
that  more  accurate  calculations  should  be  necessary  for  aerosols  with 
much  less  random  shape  and  orientation. 
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Figure  (2).  Comparison  of  non  spherical  experimental  data 
(circles)  with  figurative  Mie  phase  function  (solid  curve) 
and  two  composite  non  spherical  phase  functions  (dashed 
curves.)  The  dashed  curve  is  fitted  over  all  data,  and 
the  dashed-dotted  curve  is  fitted  only  over  data  for 
scattering  angle  greater  than  sixty  degrees.  Phase  functions 
for  shape  functions  incorporating  -.05,  -.1,  and  -.15 
were  used  in  this  fit. 
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